Introduction

To organise input data for segmented regression of age-standardised mortality rates. * Load packages * Load data * Ensure correctly labelled * Convert dates to identifiable quarters * Visualise trends

Load packages

#install.packages("pacman")
pacman::p_load(
  tidyverse,
  segmented,
  plotly
               )

Load data and assign to object “data”

data <- read_csv("Data/QuarterlyASMR1990_2018.csv")
Parsed with column specification:
cols(
  period = col_integer(),
  sex = col_character(),
  allageesprate = col_double(),
  allagestderr = col_double(),
  allagelcli = col_double(),
  allageucli = col_double(),
  under75esprate = col_double(),
  under75stderr = col_double(),
  under75lcli = col_double(),
  under75ucli = col_double(),
  totdeaths = col_integer()
)

Look at data contents

data
glimpse(data)
Observations: 342
Variables: 11
$ period         <int> 19901, 19901, 19901, 19902, 19902, 19902, 19903, 19903, 19903, 19904, 19904, 19904, 19911, 19911, 19...
$ sex            <chr> "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", "P", "F", "M", ...
$ allageesprate  <dbl> 1500.400, 2287.741, 1796.271, 1499.669, 2271.778, 1790.583, 1485.507, 2252.730, 1774.039, 1386.232, ...
$ allagestderr   <dbl> 7.788046, 13.950236, 6.854256, 7.774313, 13.869048, 6.830757, 7.730787, 13.786253, 6.788831, 7.49520...
$ allagelcli     <dbl> 1485.136, 2260.399, 1782.837, 1484.431, 2244.595, 1777.194, 1470.355, 2225.709, 1760.733, 1371.541, ...
$ allageucli     <dbl> 1515.665, 2315.084, 1809.706, 1514.907, 2298.962, 1803.971, 1500.660, 2279.751, 1787.345, 1400.923, ...
$ under75esprate <dbl> 589.2652, 1010.2307, 777.9763, 590.4938, 1002.3490, 775.1230, 588.0220, 993.3030, 770.2041, 563.1330...
$ under75stderr  <dbl> 4.994766, 7.222368, 4.242511, 4.998890, 7.183408, 4.230873, 4.985508, 7.136414, 4.212956, 4.878688, ...
$ under75lcli    <dbl> 579.4755, 996.0749, 769.6609, 580.6960, 988.2696, 766.8305, 578.2504, 979.3157, 761.9467, 553.5707, ...
$ under75ucli    <dbl> 599.0550, 1024.3866, 786.2916, 600.2917, 1016.4285, 783.4155, 597.7936, 1007.2904, 778.4615, 572.695...
$ totdeaths      <int> 34163, 31360, 65523, 34216, 31247, 65463, 33997, 31069, 65066, 31910, 29617, 61527, 31664, 29208, 60...

Separate combined dates to years and quarters

data %>%
  separate(period, into = c("year","quarter"), sep = 4) %>% 
  mutate(
    year = as.double(year),
    quarter = as.double(quarter)
  ) %>% 
  mutate(
    time = year + (quarter / 4) - 0.125
  ) %>% 
  select(time, sex, totdeaths) -> tidied_totdeaths

Visualise number of deaths by sex

tidied_totdeaths %>%
  ggplot(
    aes(
      x = time, 
      y = totdeaths,
      color = sex
    )
  ) +
  geom_line() + 
  coord_cartesian(
    ylim = c(0, 70000)
  ) -> alldeathsplot

Visualise with plotly

alldeathsplot %>%
  ggplotly()

Visualise deaths men and women only

tidied_totdeaths %>%
  filter(
    sex != "P"
  ) %>%
  ggplot(
    aes(
      x = time, 
      y = totdeaths,
      color = sex
    )
  ) +
  geom_line() + 
  coord_cartesian(
    ylim = c(0, 40000)
  )

Tidy data - years and quarters into doubles; keep only all age, assign to object “tidied_allagerates”

data %>%
  separate(period, into = c("year","quarter"), sep = 4) %>% 
  mutate(
    year = as.double(year),
    quarter = as.double(quarter)
  ) %>% 
  mutate(
    time = year + (quarter / 4) - 0.125
  ) %>% 
  select(time, sex, allageesprate, allagelcli, allageucli) -> tidied_allagerates 

Visualise ASMR with ribbon

tidied_allagerates %>%
  filter(
    sex != "P"
  ) %>%
  ggplot(
    aes(
      x = time, 
      group = sex
    )
  ) +
  geom_ribbon(
    aes(
      ymin = allagelcli,
      ymax = allageucli,
      fill = sex
      )
  ) + 
    geom_line(
    aes(y = allageesprate)
  ) + 
  coord_cartesian(
    ylim = c(0, 3000)
  )  

Segmented regression

Generate linear model for both sexes

tidied_allagerates %>% 
  rename(rate = allageesprate) %>%
  filter(sex == "M") -> allagemale
  lm(rate ~ time, data = allagemale) -> mod_male
tidied_allagerates %>%
  rename(rate = allageesprate) %>%
  filter(sex == "F") -> allagefemale
  lm(rate ~ time, data = allagefemale) -> mod_female

View linear model

summary (mod_male)

Call:
lm(formula = rate ~ time, data = allagemale)

Residuals:
   Min     1Q Median     3Q    Max 
-68.73 -38.87 -16.80  36.84 142.55 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 69951.244   1164.509   60.07   <2e-16 ***
time          -34.056      0.581  -58.62   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 51.04 on 112 degrees of freedom
Multiple R-squared:  0.9684,    Adjusted R-squared:  0.9681 
F-statistic:  3436 on 1 and 112 DF,  p-value: < 2.2e-16
summary (mod_female)

Call:
lm(formula = rate ~ time, data = allagefemale)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.519 -24.836  -3.137  17.272  89.669 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34630.7949   734.8493   47.13   <2e-16 ***
time          -16.6907     0.3666  -45.52   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 32.21 on 112 degrees of freedom
Multiple R-squared:  0.9487,    Adjusted R-squared:  0.9483 
F-statistic:  2072 on 1 and 112 DF,  p-value: < 2.2e-16

Run one break segmented regression for both sexes and view

segmented(mod_male, seg.Z = ~time, psi = 1994.125) -> maleonebreak
segmented(mod_female, seg.Z = ~time, psi = 1994.125) -> femaleonebreak
summary(maleonebreak)

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = mod_male, seg.Z = ~time, psi = 1994.125)

Estimated Break-Point(s):
     Est.   St.Err 
2013.496    0.464 

Meaningful coefficients of the linear terms:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 76334.045   1118.969  68.218   <2e-16 ***
time          -37.249      0.559 -66.637   <2e-16 ***
U1.time        39.113      5.730   6.826       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 36.76 on 110 degrees of freedom
Multiple R-Squared: 0.9839,  Adjusted R-squared: 0.9835 

Convergence attained in 2 iterations with relative change -1.957553e-16 
summary(femaleonebreak)

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = mod_female, seg.Z = ~time, psi = 1994.125)

Estimated Break-Point(s):
     Est.   St.Err 
2014.373    0.577 

Meaningful coefficients of the linear terms:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37363.1522   793.9771  47.058   <2e-16 ***
time          -18.0576     0.3966 -45.535   <2e-16 ***
U1.time        24.0687     5.4286   4.434       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 27.34 on 110 degrees of freedom
Multiple R-Squared: 0.9637,  Adjusted R-squared: 0.9627 

Convergence attained in 5 iterations with relative change 1.770497e-16 

Generate predicted values based on model, and plot alongside observed values

tidied_allagerates %>%
  filter(sex == "M") %>% 
  mutate(
    male_onebreak_pred = predict(
      maleonebreak, 
      newdata = tidied_allagerates %>% 
        filter(sex == "M") %>% 
        select(time)
      )
    ) -> malepredicted 
malepredicted %>%
  ggplot(
    aes(x = time)
  ) + 
  geom_point(
    aes(y = allageesprate)
  ) + 
  geom_line(
    aes(y = male_onebreak_pred)
  )

Merge male and female predicted values and assign to object “obs_pred_both”

Plot male and female observed and predicted

obs_pred_both %>%
  ggplot(
    aes(x = time, color = sex)
  ) + 
  geom_point(
    aes(y = observed)
  ) + 
  geom_line(
    aes(y = predicted)
  )

Generate and plot residuals based on one break model, x intercepts manually entered

obs_pred_both %>% 
  mutate(residual=observed-predicted) %>% 
  ggplot(
    aes(x = time, shape = sex, color=sex)
  ) + 
  geom_point(
    aes(y = residual)
  )+ 
  ylim(c(-100,100))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 2013.496)+
  geom_vline(xintercept = 2014.373, linetype= "dashed") 

Generate and plot residuals, plotting sexes on separate charts, x intercepts drawn direct from model

obs_pred_both %>% 
  mutate(residual=observed-predicted) %>% 
  ggplot(
    aes(x = time)
  ) + 
  geom_point(
    aes(y = residual)
  ) +
  facet_wrap(~sex) +
  ylim(c(-100,100))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = maleonebreak[["psi"]][[2]])+
    #geom_vline(xintercept = 2013.496)+
  geom_vline(xintercept = femaleonebreak[["psi"]][[2]], linetype= "dashed") 

    #geom_vline(xintercept = 2014.373, linetype= "dashed") 

Two break segmented model

segmented(mod_male, seg.Z = ~time, psi = c(1994, 2002)) -> maletwobreaks
segmented(mod_female, seg.Z = ~time, psi = c(1994, 2002)) -> femaletwobreaks
LS0tDQp0aXRsZTogIkRhdGEgbWFuYWdlbWVudCBmb3Igc2VnbWVudGVkIHJlZ3Jlc3Npb24iDQpvdXRwdXQ6DQogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQNCiAgaHRtbF9kb2N1bWVudDoNCiAgICBkZl9wcmludDogcGFnZWQNCiAgd29yZF9kb2N1bWVudDogZGVmYXVsdA0KLS0tDQojIEludHJvZHVjdGlvbg0KDQpUbyBvcmdhbmlzZSBpbnB1dCBkYXRhIGZvciBzZWdtZW50ZWQgcmVncmVzc2lvbiBvZiBhZ2Utc3RhbmRhcmRpc2VkIG1vcnRhbGl0eSByYXRlcy4gDQoqIExvYWQgcGFja2FnZXMNCiogTG9hZCBkYXRhDQoqIEVuc3VyZSBjb3JyZWN0bHkgbGFiZWxsZWQNCiogQ29udmVydCBkYXRlcyB0byBpZGVudGlmaWFibGUgcXVhcnRlcnMNCiogVmlzdWFsaXNlIHRyZW5kcw0KDQojIyBMb2FkIHBhY2thZ2VzDQoNCmBgYHtyfQ0KI2luc3RhbGwucGFja2FnZXMoInBhY21hbiIpDQpwYWNtYW46OnBfbG9hZCgNCiAgdGlkeXZlcnNlLA0KICBzZWdtZW50ZWQsDQogIHBsb3RseQ0KICAgICAgICAgICAgICAgKQ0KYGBgDQoNCiMjIExvYWQgZGF0YSBhbmQgYXNzaWduIHRvIG9iamVjdCAiZGF0YSINCg0KYGBge3J9DQpkYXRhIDwtIHJlYWRfY3N2KCJEYXRhL1F1YXJ0ZXJseUFTTVIxOTkwXzIwMTguY3N2IikNCmBgYA0KDQojIyBMb29rIGF0IGRhdGEgY29udGVudHMNCg0KYGBge3J9DQpkYXRhDQpnbGltcHNlKGRhdGEpDQpgYGANCg0KIyMgU2VwYXJhdGUgY29tYmluZWQgZGF0ZXMgdG8geWVhcnMgYW5kIHF1YXJ0ZXJzDQoNCmBgYHtyfQ0KZGF0YSAlPiUNCiAgc2VwYXJhdGUocGVyaW9kLCBpbnRvID0gYygieWVhciIsInF1YXJ0ZXIiKSwgc2VwID0gNCkgJT4lIA0KICBtdXRhdGUoDQogICAgeWVhciA9IGFzLmRvdWJsZSh5ZWFyKSwNCiAgICBxdWFydGVyID0gYXMuZG91YmxlKHF1YXJ0ZXIpDQogICkgJT4lIA0KICBtdXRhdGUoDQogICAgdGltZSA9IHllYXIgKyAocXVhcnRlciAvIDQpIC0gMC4xMjUNCiAgKSAlPiUgDQogIHNlbGVjdCh0aW1lLCBzZXgsIHRvdGRlYXRocykgLT4gdGlkaWVkX3RvdGRlYXRocw0KYGBgDQoNCg0KIyMgVmlzdWFsaXNlIG51bWJlciBvZiBkZWF0aHMgYnkgc2V4DQoNCmBgYHtyfQ0KdGlkaWVkX3RvdGRlYXRocyAlPiUNCiAgZ2dwbG90KA0KICAgIGFlcygNCiAgICAgIHggPSB0aW1lLCANCiAgICAgIHkgPSB0b3RkZWF0aHMsDQogICAgICBjb2xvciA9IHNleA0KICAgICkNCiAgKSArDQogIGdlb21fbGluZSgpICsgDQogIGNvb3JkX2NhcnRlc2lhbigNCiAgICB5bGltID0gYygwLCA3MDAwMCkNCiAgKSAtPiBhbGxkZWF0aHNwbG90DQoNCmBgYA0KDQojIyBWaXN1YWxpc2Ugd2l0aCBwbG90bHkNCg0KYGBge3J9DQphbGxkZWF0aHNwbG90ICU+JQ0KICBnZ3Bsb3RseSgpDQoNCmBgYA0KIyMgVmlzdWFsaXNlIGRlYXRocyBtZW4gYW5kIHdvbWVuIG9ubHkNCg0KYGBge3J9DQp0aWRpZWRfdG90ZGVhdGhzICU+JQ0KICBmaWx0ZXIoDQogICAgc2V4ICE9ICJQIg0KICApICU+JQ0KICBnZ3Bsb3QoDQogICAgYWVzKA0KICAgICAgeCA9IHRpbWUsIA0KICAgICAgeSA9IHRvdGRlYXRocywNCiAgICAgIGNvbG9yID0gc2V4DQogICAgKQ0KICApICsNCiAgZ2VvbV9saW5lKCkgKyANCiAgY29vcmRfY2FydGVzaWFuKA0KICAgIHlsaW0gPSBjKDAsIDQwMDAwKQ0KICApDQoNCmBgYA0KIyMgVGlkeSBkYXRhIC0geWVhcnMgYW5kIHF1YXJ0ZXJzIGludG8gZG91Ymxlczsga2VlcCBvbmx5IGFsbCBhZ2UsIGFzc2lnbiB0byBvYmplY3QgInRpZGllZF9hbGxhZ2VyYXRlcyINCg0KYGBge3J9DQpkYXRhICU+JQ0KICBzZXBhcmF0ZShwZXJpb2QsIGludG8gPSBjKCJ5ZWFyIiwicXVhcnRlciIpLCBzZXAgPSA0KSAlPiUgDQogIG11dGF0ZSgNCiAgICB5ZWFyID0gYXMuZG91YmxlKHllYXIpLA0KICAgIHF1YXJ0ZXIgPSBhcy5kb3VibGUocXVhcnRlcikNCiAgKSAlPiUgDQogIG11dGF0ZSgNCiAgICB0aW1lID0geWVhciArIChxdWFydGVyIC8gNCkgLSAwLjEyNQ0KICApICU+JSANCiAgc2VsZWN0KHRpbWUsIHNleCwgYWxsYWdlZXNwcmF0ZSwgYWxsYWdlbGNsaSwgYWxsYWdldWNsaSkgLT4gdGlkaWVkX2FsbGFnZXJhdGVzIA0KYGBgDQoNCiMjVmlzdWFsaXNlIEFTTVIgd2l0aCByaWJib24NCg0KYGBge3J9DQp0aWRpZWRfYWxsYWdlcmF0ZXMgJT4lDQogIGZpbHRlcigNCiAgICBzZXggIT0gIlAiDQogICkgJT4lDQogIGdncGxvdCgNCiAgICBhZXMoDQogICAgICB4ID0gdGltZSwgDQogICAgICBncm91cCA9IHNleA0KICAgICkNCiAgKSArDQoNCiAgZ2VvbV9yaWJib24oDQogICAgYWVzKA0KICAgICAgeW1pbiA9IGFsbGFnZWxjbGksDQogICAgICB5bWF4ID0gYWxsYWdldWNsaSwNCiAgICAgIGZpbGwgPSBzZXgNCiAgICAgICkNCiAgKSArIA0KICAgIGdlb21fbGluZSgNCiAgICBhZXMoeSA9IGFsbGFnZWVzcHJhdGUpDQogICkgKyANCiAgY29vcmRfY2FydGVzaWFuKA0KICAgIHlsaW0gPSBjKDAsIDMwMDApDQogICkgIA0KYGBgDQoNCiMgU2VnbWVudGVkIHJlZ3Jlc3Npb24gDQojI0dlbmVyYXRlIGxpbmVhciBtb2RlbCBmb3IgYm90aCBzZXhlcw0KDQpgYGB7cn0NCnRpZGllZF9hbGxhZ2VyYXRlcyAlPiUgDQogIHJlbmFtZShyYXRlID0gYWxsYWdlZXNwcmF0ZSkgJT4lDQogIGZpbHRlcihzZXggPT0gIk0iKSAtPiBhbGxhZ2VtYWxlDQogIGxtKHJhdGUgfiB0aW1lLCBkYXRhID0gYWxsYWdlbWFsZSkgLT4gbW9kX21hbGUNCg0KdGlkaWVkX2FsbGFnZXJhdGVzICU+JQ0KICByZW5hbWUocmF0ZSA9IGFsbGFnZWVzcHJhdGUpICU+JQ0KICBmaWx0ZXIoc2V4ID09ICJGIikgLT4gYWxsYWdlZmVtYWxlDQogIGxtKHJhdGUgfiB0aW1lLCBkYXRhID0gYWxsYWdlZmVtYWxlKSAtPiBtb2RfZmVtYWxlDQoNCmBgYA0KDQojI1ZpZXcgbGluZWFyIG1vZGVsIA0KDQpgYGB7cn0NCnN1bW1hcnkgKG1vZF9tYWxlKQ0Kc3VtbWFyeSAobW9kX2ZlbWFsZSkNCg0KYGBgDQoNCiMjUnVuIG9uZSBicmVhayBzZWdtZW50ZWQgcmVncmVzc2lvbiBmb3IgYm90aCBzZXhlcyBhbmQgdmlldw0KDQpgYGB7cn0NCnNlZ21lbnRlZChtb2RfbWFsZSwgc2VnLlogPSB+dGltZSwgcHNpID0gMTk5NC4xMjUpIC0+IG1hbGVvbmVicmVhaw0Kc2VnbWVudGVkKG1vZF9mZW1hbGUsIHNlZy5aID0gfnRpbWUsIHBzaSA9IDE5OTQuMTI1KSAtPiBmZW1hbGVvbmVicmVhaw0Kc3VtbWFyeShtYWxlb25lYnJlYWspDQpzdW1tYXJ5KGZlbWFsZW9uZWJyZWFrKQ0KDQpgYGANCg0KIyNHZW5lcmF0ZSBwcmVkaWN0ZWQgdmFsdWVzIGJhc2VkIG9uIG1vZGVsLCBhbmQgcGxvdCBhbG9uZ3NpZGUgb2JzZXJ2ZWQgdmFsdWVzDQoNCmBgYHtyfQ0KdGlkaWVkX2FsbGFnZXJhdGVzICU+JQ0KICBmaWx0ZXIoc2V4ID09ICJNIikgJT4lIA0KICBtdXRhdGUoDQogICAgbWFsZV9vbmVicmVha19wcmVkID0gcHJlZGljdCgNCiAgICAgIG1hbGVvbmVicmVhaywgDQogICAgICBuZXdkYXRhID0gdGlkaWVkX2FsbGFnZXJhdGVzICU+JSANCiAgICAgICAgZmlsdGVyKHNleCA9PSAiTSIpICU+JSANCiAgICAgICAgc2VsZWN0KHRpbWUpDQogICAgICApDQogICAgKSAtPiBtYWxlcHJlZGljdGVkIA0KbWFsZXByZWRpY3RlZCAlPiUNCiAgZ2dwbG90KA0KICAgIGFlcyh4ID0gdGltZSkNCiAgKSArIA0KICBnZW9tX3BvaW50KA0KICAgIGFlcyh5ID0gYWxsYWdlZXNwcmF0ZSkNCiAgKSArIA0KICBnZW9tX2xpbmUoDQogICAgYWVzKHkgPSBtYWxlX29uZWJyZWFrX3ByZWQpDQogICkNCmBgYA0KDQojI01lcmdlIG1hbGUgYW5kIGZlbWFsZSBwcmVkaWN0ZWQgdmFsdWVzIGFuZCBhc3NpZ24gdG8gb2JqZWN0ICJvYnNfcHJlZF9ib3RoIg0KDQpgYGB7cn0NCmZlbWFsZXByZWRpY3RlZCAlPiUgDQogIHNlbGVjdCgNCiAgICB0aW1lLCBzZXgsIG9ic2VydmVkID0gYWxsYWdlZXNwcmF0ZSwNCiAgICBsb3dlciA9IGFsbGFnZWxjbGksIHVwcGVyID0gYWxsYWdldWNsaSwNCiAgICBwcmVkaWN0ZWQgPSBmZW1hbGVfb25lYnJlYWtfcHJlZA0KICApICU+JSANCiAgYmluZF9yb3dzKA0KICAgIG1hbGVwcmVkaWN0ZWQgJT4lIA0KICAgICAgc2VsZWN0KA0KICAgICAgICB0aW1lLCBzZXgsIG9ic2VydmVkID0gYWxsYWdlZXNwcmF0ZSwNCiAgICAgICAgbG93ZXIgPSBhbGxhZ2VsY2xpLCB1cHBlciA9IGFsbGFnZXVjbGksDQogICAgICAgIHByZWRpY3RlZCA9IG1hbGVfb25lYnJlYWtfcHJlZA0KICAgICAgKSANCiAgKSAtPiBvYnNfcHJlZF9ib3RoDQoNCg0KYGBgDQoNCiMjUGxvdCBtYWxlIGFuZCBmZW1hbGUgb2JzZXJ2ZWQgYW5kIHByZWRpY3RlZA0KDQpgYGB7cn0NCm9ic19wcmVkX2JvdGggJT4lDQogIGdncGxvdCgNCiAgICBhZXMoeCA9IHRpbWUsIGNvbG9yID0gc2V4KQ0KICApICsgDQogIGdlb21fcG9pbnQoDQogICAgYWVzKHkgPSBvYnNlcnZlZCkNCiAgKSArIA0KICBnZW9tX2xpbmUoDQogICAgYWVzKHkgPSBwcmVkaWN0ZWQpDQogICkNCmBgYA0KDQojI0dlbmVyYXRlIGFuZCBwbG90IHJlc2lkdWFscyBiYXNlZCBvbiBvbmUgYnJlYWsgbW9kZWwsIHggaW50ZXJjZXB0cyBtYW51YWxseSBlbnRlcmVkDQoNCmBgYHtyfQ0Kb2JzX3ByZWRfYm90aCAlPiUgDQogIG11dGF0ZShyZXNpZHVhbD1vYnNlcnZlZC1wcmVkaWN0ZWQpICU+JSANCiAgZ2dwbG90KA0KICAgIGFlcyh4ID0gdGltZSwgc2hhcGUgPSBzZXgsIGNvbG9yPXNleCkNCiAgKSArIA0KICBnZW9tX3BvaW50KA0KICAgIGFlcyh5ID0gcmVzaWR1YWwpDQogICkrIA0KICB5bGltKGMoLTEwMCwxMDApKSsNCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMCkrDQogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDIwMTMuNDk2KSsNCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMjAxNC4zNzMsIGxpbmV0eXBlPSAiZGFzaGVkIikgDQoNCmBgYA0KDQojI0dlbmVyYXRlIGFuZCBwbG90IHJlc2lkdWFscywgcGxvdHRpbmcgc2V4ZXMgb24gc2VwYXJhdGUgY2hhcnRzLCB4IGludGVyY2VwdHMgZHJhd24gZGlyZWN0IGZyb20gbW9kZWwNCg0KYGBge3J9DQpvYnNfcHJlZF9ib3RoICU+JSANCiAgbXV0YXRlKHJlc2lkdWFsPW9ic2VydmVkLXByZWRpY3RlZCkgJT4lIA0KICBnZ3Bsb3QoDQogICAgYWVzKHggPSB0aW1lKQ0KICApICsgDQogIGdlb21fcG9pbnQoDQogICAgYWVzKHkgPSByZXNpZHVhbCkNCiAgKSArDQogIGZhY2V0X3dyYXAofnNleCkgKw0KICB5bGltKGMoLTEwMCwxMDApKSsNCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMCkrDQogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IG1hbGVvbmVicmVha1tbInBzaSJdXVtbMl1dKSsNCiAgICAjZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMjAxMy40OTYpKw0KICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSBmZW1hbGVvbmVicmVha1tbInBzaSJdXVtbMl1dLCBsaW5ldHlwZT0gImRhc2hlZCIpIA0KICAgICNnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAyMDE0LjM3MywgbGluZXR5cGU9ICJkYXNoZWQiKSANCg0KDQpgYGANCg0KIyNUd28gYnJlYWsgc2VnbWVudGVkIG1vZGVsICAgIA0KDQpgYGB7cn0NCnNlZ21lbnRlZChtb2RfbWFsZSwgc2VnLlogPSB+dGltZSwgcHNpID0gYygxOTk0LCAyMDAyKSkgLT4gbWFsZXR3b2JyZWFrcw0Kc2VnbWVudGVkKG1vZF9mZW1hbGUsIHNlZy5aID0gfnRpbWUsIHBzaSA9IGMoMTk5NCwgMjAwMikpIC0+IGZlbWFsZXR3b2JyZWFrcw0KDQpgYGA=